Understanding patterns of species abundance and persistence in IBC-grass

1 Monoculture biomass vs. traits

IBC_grass.monocultures <- IBC_grass %>%
    filter(Ninitial == 1,
           Year == 200) %>%
    rename(monoculture.biomass = Biomass) %>%
    select(-Ninitial, -Stage, -Year, -Rep) %>%
    pivot_longer(c(-SpeciesID, -monoculture.biomass),
                 names_to = "trait", 
                 values_to = "trait.value")

ggplot(IBC_grass.monocultures) +
    geom_boxplot(aes(x = as.factor(trait.value),
                     y = monoculture.biomass)) +
    facet_wrap(facets = vars(trait),
               ncol = 4,
               scales = "free_x") +
    labs(x = "Trait value",
         y = "Monoculture biomass") +
    theme_bw()

IBC_grass.monocultures <- IBC_grass %>%
    filter(Ninitial == 1,
           Year %in% seq(180, 200)) %>%
    rename(monoculture.biomass = Biomass) %>%
    mutate(id = row_number()) %>%
    mutate_if(is.character, as.factor) %>%
    select(-SpeciesID, -Ninitial, -Stage, -Rep) 

train <- IBC_grass.monocultures %>% sample_frac(.70)
test <- anti_join(IBC_grass.monocultures, train, by = 'id')

train <- train %>% select(-id)
test <- test %>% select(-id)

rf <- randomForest(data = as.data.frame(train),
                   monoculture.biomass ~ .,
                   importance = TRUE)

pred <- data.frame(pred = predict(rf, test))
title = paste("correlation: ", round(cor(pred, test$monoculture.biomass)[[1]], 2),
              "  |  mean squared error: ", round(mean(rf$mse), 2),
              "  |  R-squared: ", round(mean(rf$rsq), 2),
              sep = "")

varImpPlot(rf, main = title)

IBC_grass.monocultures <- IBC_grass %>%
    filter(Ninitial == 1,
           Year == 200) %>%
    select(-Ninitial, -Stage, -Year, -Rep)

ggplot(IBC_grass.monocultures) +
    geom_point(aes(x = Gmax, 
                   y = MaxMass,
                   color = Biomass),
               size = 4) +
    scale_color_viridis() +
    labs(color = "Monoculture\nbiomass") +
    theme_bw()

IBC_grass.monocultures <- IBC_grass %>%
    filter(Ninitial == 1,
           Year == 200) %>%
    rename(monoculture.biomass = Biomass) %>%
    select(-Ninitial, -Stage, -Year, -Rep) 

IBC_grass.32 <- IBC_grass %>%
    filter(Ninitial == 32,
           Year == 200) %>%
    select(SpeciesID, Biomass)

IBC_grass.combined <- inner_join(IBC_grass.monocultures, IBC_grass.32)

ggplot(IBC_grass.combined) +
    geom_point(aes(x = monoculture.biomass,
                   y = Biomass,
                   color = SpeciesID)) +
    labs(x = "Monoculture biomass",
         y = "Biomass in 32-species mixture") +
    theme_bw() + 
    theme(legend.position = "none")

IBC_grass.32 <- IBC_grass %>%
    filter(Ninitial == 32,
           Year == 200) %>%
    select(SpeciesID, Biomass, LMR, MaxMass, Gmax, SLA, meanSpacerLength)

IBC_grass.32 <- IBC_grass.32 %>%
    pivot_longer(c(-SpeciesID, -Biomass),
                 names_to = "trait",
                 values_to = "trait.value")

ggplot(IBC_grass.32) +
    geom_point(aes(x = trait.value,
                   y = Biomass,
                   color = SpeciesID)) +
    facet_wrap(facets = vars(trait),
               ncol = 4,
               scales = "free_x") +
    labs(x = "Trait value",
         y = "Biomass in 32-species mixture") +
    theme_bw() +
    theme(legend.position = "none")

IBC_grass.32 <- IBC_grass %>%
    filter(Ninitial == 32,
           Year %in% seq(180, 200))

IBC_grass.32 <- IBC_grass.32 %>%
    ungroup() %>%
    mutate(id = row_number()) %>%
    mutate_if(is.character, as.factor) %>%
    select(-Ninitial, -Rep, -Stage, -SpeciesID) 

train <- IBC_grass.32 %>% sample_frac(.70)
test <- anti_join(IBC_grass.32, train, by = 'id')

train <- train %>% select(-id)
test <- test %>% select(-id)

rf <- randomForest(data = as.data.frame(train),
                   Biomass ~ .,
                   importance = TRUE)

pred <- data.frame(pred = predict(rf, test))

title = paste("correlation: ", round(cor(pred, test$Biomass)[[1]], 2),
              "  |  mean squared error: ", round(mean(rf$mse), 2),
              "  |  R-squared: ", round(mean(rf$rsq), 2),
              sep = "")

varImpPlot(rf, main = title)

IBC_grass.32 <- IBC_grass %>%
    filter(Ninitial == 32,
           Year == 200) %>%
    group_by(SpeciesID, LMR, MaxMass, Gmax, SLA, meanSpacerLength) %>%
    summarise(Biomass = mean(Biomass)) %>%
    filter(Biomass > 0)

plot_ly(IBC_grass.32,
        x = ~Gmax, 
        y = ~MaxMass,
        z = ~SLA,
        color = ~log(Biomass),
        mode = "markers",
        type = "scatter3d")

IBC_grass.monocultures <- IBC_grass %>%
    filter(Ninitial == 1,
           Year == 200) %>%
    rename(monoculture.biomass = Biomass) %>%
    select(-Ninitial, -Stage, -Year, -Rep) 

IBC_grass.32 <- IBC_grass %>%
    filter(Ninitial == 32)

IBC_grass.32 <- IBC_grass.32 %>%
    group_by(SpeciesID, Rep) %>%
    arrange(Year) %>%
    mutate(isDead = (Biomass == 0)) %>%
    filter(isDead | (!isDead & Year == max(Year))) %>%
    filter(Year == first(Year))

IBC_grass.32 <- IBC_grass.32 %>%
    group_by(SpeciesID) %>%
    summarise(Year = mean(Year))

IBC_grass.combined <- inner_join(IBC_grass.monocultures, IBC_grass.32)

ggplot(IBC_grass.combined) +
    geom_point(aes(x = monoculture.biomass,
                   y = Year,
                   color = SpeciesID)) +
    labs(x = "Monoculture biomass",
         y = "Extinction year") +
    theme_bw() +
    theme(legend.position = "none")

IBC_grass.32 <- IBC_grass %>%
    filter(Ninitial == 32)

IBC_grass.32 <- IBC_grass.32 %>%
    group_by(SpeciesID, Rep) %>%
    arrange(Year) %>%
    mutate(isDead = (Biomass == 0)) %>%
    filter(isDead | (!isDead & Year == max(Year))) %>%
    filter(Year == first(Year))

IBC_grass.32 <- IBC_grass.32 %>%
    group_by(SpeciesID, LMR, MaxMass, Gmax, SLA, meanSpacerLength) %>%
    summarise(Year = mean(Year))

IBC_grass.32 <- IBC_grass.32 %>%
    pivot_longer(c(-SpeciesID, -Year),
                 names_to = "trait",
                 values_to = "trait.value")

ggplot(IBC_grass.32) +
    geom_point(aes(x = trait.value,
                   y = Year,
                   color = SpeciesID)) +
    facet_wrap(facets = vars(trait),
               ncol = 4,
               scales = "free") +
    labs(x = "Trait value",
         y = "Time to extinction") +
    theme_bw() +
    theme(legend.position = "none")

IBC_grass.32 <- IBC_grass %>%
    filter(Ninitial == 32)

IBC_grass.32 <- IBC_grass.32 %>%
    group_by(SpeciesID, Rep) %>%
    arrange(Year) %>%
    mutate(isDead = (Biomass == 0)) %>%
    filter(isDead | (!isDead & Year == max(Year))) %>%
    filter(Year == first(Year)) %>%
    rename(extinction.year = Year)

IBC_grass.32 <- IBC_grass.32 %>%
    ungroup() %>%
    mutate(id = row_number()) %>%
    mutate_if(is.character, as.factor) %>%
    select(-SpeciesID, -Ninitial, -Stage, -Rep, -Biomass, -isDead) 

train <- IBC_grass.32 %>% sample_frac(.70)
test <- anti_join(IBC_grass.32, train, by = 'id')

train <- train %>% select(-id)
test <- test %>% select(-id)

rf <- randomForest(data = as.data.frame(train),
                   extinction.year ~ .,
                   importance = TRUE)

pred <- data.frame(pred = predict(rf, test))

title = paste("correlation: ", round(cor(pred, test$extinction.year)[[1]], 2),
              "  |  mean squared error: ", round(mean(rf$mse), 2),
              "  |  R-squared: ", round(mean(rf$rsq), 2),
              sep = "")

varImpPlot(rf, main = title)

IBC_grass.32 <- IBC_grass %>%
    filter(Ninitial == 32)

IBC_grass.32 <- IBC_grass.32 %>%
    group_by(SpeciesID, Rep) %>%
    arrange(Year) %>%
    mutate(isDead = (Biomass == 0)) %>%
    filter(isDead | (!isDead & Year == max(Year))) %>%
    filter(Year == first(Year))

IBC_grass.32 <- IBC_grass.32 %>%
    group_by(SpeciesID, LMR, MaxMass, Gmax, SLA, meanSpacerLength) %>%
    summarise(Year = mean(Year))

ggplot(IBC_grass.32) +
    geom_point(aes(x = Gmax,
                   y = MaxMass,
                   color = Year)) +
    labs(color = "Extinction\nyear") +
    scale_color_viridis() +
    theme_bw()

IBC_grass.monocultures <- IBC_grass %>%
    filter(Ninitial == 1,
           Year %in% c(100, 200)) %>%
    rename(monoculture.biomass = Biomass) %>%
    select(-Ninitial, -Year, -Rep) 
    
IBC_grass.mixture <- IBC_grass %>%
    filter(Ninitial == 64,
           Year %in% c(100, 200)) %>%
    rename(mixture.biomass = Biomass) %>%
    select(-Ninitial, -Year, -Rep) 

IBC_grass.combined <- inner_join(IBC_grass.monocultures, IBC_grass.mixture)

ggplot(IBC_grass.combined) +
    geom_point(aes(x = monoculture.biomass,
                   y = mixture.biomass,
                   color = SpeciesID)) +
    facet_grid(cols = vars(Stage)) +
    labs(x = "Monoculture biomass",
         y = "Mixture biomass") + 
    theme_bw() +
    theme(legend.position = "none")

IBC_grass.monocultures <- IBC_grass %>%
    filter(Ninitial == 1,
           Year == 200) %>%
    rename(monoculture.biomass = Biomass) %>%
    select(-Ninitial, -Year, -Rep) 
    
IBC_grass.64 <- IBC_grass %>%
    filter(Ninitial == 64) %>%
    select(-Ninitial, -Rep)

IBC_grass.64 <- IBC_grass.64 %>%
    group_by(SpeciesID) %>%
    arrange(Year) %>%
    mutate(isDead = (Biomass == 0)) %>%
    filter(isDead | (!isDead & Year == max(Year))) %>%
    filter(Year == first(Year))
    
IBC_grass.combined <- inner_join(IBC_grass.monocultures, IBC_grass.64)

ggplot(IBC_grass.combined) +
    geom_point(aes(x = monoculture.biomass,
                   y = Year,
                   color = SpeciesID)) +
    labs(x = "Monoculture biomass",
         y = "Extinction year") +
    theme_bw() +
    theme(legend.position = "none")